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ABSTRACT 

The  computer  program  CHIEF  'was  designed  bySchenckand  Barach 
of  the  Naval  Undersea  Research  and  Development  Center  to  obtain  ap¬ 
proximate  solutions  to  exterior  steady-state  acoustic  radiation  problems 
for  surfaces  of  arbitrary  shape  vibrating  with  a  prescribed  normal 
velocity.  To  test  its  capabilities  as  a  research  tool,  CHIEF  was  applied 
to  several  problems  for  which  accurate  answers  have  been  obtained 
using  harmonic  expansions.  The  accuracy  and  computation  time  of  the 
results  using  CHIEF  are  discussed  in  terms  of  the  surface  subdivision 
scheme  and  the  number  of  Gaussian  quadrature  points  used  to  evaluate 
the  Helmholtz  integrals.  In  addition,  new  input  subroutines  to  CHIEF 
which  provide  numen.us  geometrical  options  are  also  discussed. 


PROBLEM  STATUS 

This  is  an  interim  report  on  a  continuing  NRL  Problem, 
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A  TEST  OF  THE  CAPABILITIES  OF  CHIEF 
IN  THE  NUMERICAL  CALCULATION  OF  ACOUSTIC 
RADIATION  FROM  ARBITRARY  SURFACES 


INTRODUCTION 

The  computer  program  CHIEF  as  developed  by  Schenck  (1)  and  Barach  (?-)  is 
designed  to  obtain  approximate  solutions  to  exterior  steady- state  acoiistic  radiation 
problems  for  surfaces  of  arbitrary  shape  vibrating  with  a  prescribed  normal  velocity. 

To  test  its  capabilities  as  a  research  tool,  CHIFF  was  applied  to  several  problems  for 
which  accurate  answers  have  been  obtained  using  harmonic  expansions. 

The  first  examples  involve  the  acoustic  radiation  from  an  oblate  spheroid  whose  top 
half  is  vibrating  with  unit  normal  velocity  and  whose  bottom  half  is  rigid.  The  acoustic 
radiation  impedance  of  the  top  half  and  the  mutual  radiation  impedance  coefficient  between 
the  two  halves  of  the  spheroid  obtained  using  CHIEF  are  compared  with  results  obtained 
using  a  harmonic  expansion  in  oblate  spheroidal  harmonic  functions.  A  Fortran  com¬ 
puter  program  called  OBRAD  (3)  was  used  to  accurately  evaluate  the  necessary  sphe¬ 
roidal  functions. 

An  example  designed  to  test  CHIEF'S  ability  to  handle  a  multiple  surface  is  that  of 
the  radiation  of  a  uniformly  pulsating  sphere  in  the  presence  of  a  similar  stationary 
sphere.  Here  the  results  from  CHIEF  are  compared  with  those  of  New  (4),  who  has 
obtained  accurate  values  for  both  the  near-field  and  the  far- field  pressures  using  a 
harmonic  expansion  in  terms  of  spherical  functions. 

The  accuracy  and  computation  time  of  the  results  using  CHIEF  are  discussed  in 
terms  of  the  surface  subdivision  scheme  and  the  number  of  Gaussian  quadrature  points 
used  to  evaluate  the  non-self  Helmholtz  integrals.  In  addition,  new  input  subroutines  to 
CHIEF  which  provide  numerous  geometrical  options  are  discussed.  A  computer  print¬ 
out  of  these  subroutines  is  given  in  the  appendix. 


REVIEW  OF  THE  COMBINED  HELMHOLTZ  INTEGRAL 
EQUATION  FORMULATION  (CHIEF) 

Consider  a  finite  region  bounded  by  the  regular,  closed  surface  S,  as  shown  in 
Fig.  1.  Lei  an  arbitrary  point  on  the  surface  be  denoted  by  e.  The  region  exterior  to 
S  is  assumed  to  be  filled  with  an  ideal,  homogeneous  fluid  of  density ..  and  sound  speed 
c .  Let  an  arbitrary  point  in  this  exterior  region  be  denoted  by  x .  The  surface  S  is 
vibrating  at  an  angular  frequency  w  with  a  known  normal  velocity  distribution  v{Q)  .  The 
steady-state  pressure  p(x)  may  be  obtained  by  solving  the  Helmholtz  scalar  ^^ave  equa¬ 
tion, 


(v2+  k^)  p  (X)  =  0, 


(1) 


where  k  r  o^'c  .  The  time  dependence 


c'"*  has  been  suppressed. 
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Fig.  1.  -  An  arbitrary  surface  S 


The  solution  to  Eq.  (1)  must  be  finite  and  must  satisfy  the  radiation  condition  at 
infinity, 


where  Sr  is  the  surface  of  a  sphere  of  radius  R  surrounding  the  surface  s.  In  addition, 
the  pressure  p(i)  must  satisfy  the  boundary  condition  on  s, 

Sp 

~=-icupv(§),  (3) 

where  denotes  the  outward  normal  derivative  evaluated  at  the  surface  point  i.  The 
solution  to  Eq.  (1)  which  satisfies  the  boundary  conditions,  Eq.  (2)  and  Eq.  (3),  is  given 
by 

-  J  f  -  a  ?)1  p-ikdCx.f) 

P(x)  =  J  P(§)  nKiriTT 

where  d(i,  |)  is  the  distance  between  the  exterior  point  x  and  the  surface  point  I. 

This  expression  allows  p(x)  to  be  evaluated  when  the  pressure  on  the  surface  p(|) 
is  known.  In  order  to  obtain  p(C)t  the  point  of  observation  5  is  allowed  to  approach  the 
surface.  When  the  limiting  process  is  properly  performed,  one  obtains  the  surface 
Helmholtz  integral  formulation. 


p(|-,  =  i,  f  p(5)  =  as<i) 

2nJ  ^  <1(5',!)  J  d(r,  §) 


If  the  field  point  is  taken  interior  to  the  surface  s,  one  obtains  the  interior  Helmholtz 
integral  formulation, 

)  r  -  3  Fp-iMy  .5)1  .  I) 

" "  4;-  J  '■«>  "d<or  •  <“> 


where  y  is  the  interior  point. 
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The  CHIEF  program  solves  for  the  surface  pressure  p{|)  using  Eq.  (5)  and,  if 
necessary,  Eq.  (6). 

In  order  to  solve  Eq.  (5),  the  surface  is  subdivided  into  small  areas.  The  normal 
particle  velocity  is  chosen  to  be  constant  over  each  subdivision,  and  it  is  assumed  that 
the  pressure  io  constant  over  each  subdivision.  The  latter  assumption  is  an  approxima¬ 
tion  that  will  be  good  only  if  the  true  pressure  does  not  vary  much  over  each  subdivision. 
The  approximation  can  be  improved,  if  necessary,  by  further  subdivision  of  the  surface. 

The  surface  integrals  in  Eq.  (5)  can  now  be  broken  up  into  integrations  over  each 
subdivision  ,  giving 


(7) 


where  and  are  the  pressure  and  normal  particle  velocity  for  the  subdivision  s^j . 

If  the  observation  point  |'  is  chosen  to  be  on  S„,  one  obtains  a  set  of  simultaneous 
equations  in  the  unknown  pressures. 


=  (8) 
/3  fi 


where  and  are  given  by 


Kp  - 


0/3 


j  cn^  L  d(5„.  J 


dS(l^) , 


(9) 


where  is  the  Kronecker  delta, 
and 


'5^/3=  / 


ip) 


(10) 


The  Helmholtz  integrals  given  in  Eqs.  (9)  and  (10)  are  numerically  evaluated  in 
CHIEF  using  a  Gaussian  quadrature  over  both  surface  coordinates.  For  the  so  called 
non-self  Helmholtz  integrals,  when  a  /  /?,  the  user  of  CHIEF  must  input  the  number  of 
quadrature  points  used  to  evaluate  the  integral.  The  self  integrals,  i.e. ,  when  a  =  /?, 
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are  automatically  evaluated  by  first  subdividing  the  integrals  into  four  pieces  and  then 
evaluating  each  piece  using  16-point  Gaussian  quadratures  over  both  surface  coordinates. 
Obviously,  the  choice  of  the  number  of  quadrature  points  to  be  used  for  the  non-stlf 
integrals  is  influenced  by  the  surface  subdivision  scheme. 

If  the  frequency  is  near  one  of  the  characteristic  frequencies  for  the  interior  homo¬ 
geneous  Dirichlet  problem,  the  simultaneous  equations,  Eq,  (8),  will  not  yield  the  cor- 
re^  ’  surface  pre.csuros.  In  this  case  advUtional  equations  are  obtained  using  the  interior 
Helmholtz  integral  formulation  for  carciully  chosen  interior  points.  These  equations 
add  no  more  unknowns  but  provide  aud*'  ional  equations  of  constraint  for  the  surface 
pressures  leading  to  an  overdetermined  system.  Since  only  the  correct  set  of  surface 
pressures  satisfies  these  additional  equations,  their  addition  tends  to  force  the  solution 
to  the  desired  one.  The  clioice  of  the  number  and  the  location  of  the  interior  points  is  a 
difficult  problem.  The  examples  discussed  in  this  paper’will  not  require  any  interior 
points.  This  will  allow  an  uncluttered  examination  of  the  more  basic  question  confronting 
the  user  of  CHIEF:  how  do  both  the  accuracy  and  the  computation  time  depend  on  the 
surface  subdivision  scheme  and  the  number  of  quadrature  points  ? 


COMPARISON  OF  RESULTS 

Example  I.  Oblate  Spheroid 

The  first  example  is  that  of  radiation  from  an  oblate  spheroid,  as  shown  in  Fig.  2. 
Here  the  radial  coordinate  specifying  the  oblate  spheroidal  surface,  is  chosen  equal 
to  0.2.  The  ratio  of  major  to  minor  axes  is  then  given  by  v/(f^^TT)/^=  5.0.  The  param¬ 
eter  l»,  defined  to  be  v  times  the  ratio  of  the  distance  between  the  foci  of  the  elliptical 
cross  section  to  the  wavelength,  is  a  measure  of  the  acoustical  size.  Here  h  is  chosen 
equal  to  unity  to  insure  that  the  frequency  is  well  below  the  lowest  characteristic  fre¬ 
quency  for  the  interior  homogeneous  Dirichlet  problem.  For  this  example  the  top  half 
of  the  surface  is  specified  to  be  vibrating  with  unit  normal  velocity.  Of  interest  are  the 
self  acoustic  r'ldiation  impedance  coefficient  of  the  top  half  and  the  mutual  acoustic 
radiation  impedance  coefficient  between  the  two  halves.  These  were  first  calculated 
analytically  using  a  harmonic  expansion  in  oblate  spheroidal  wave  functions.  Sufficient 
spheroidal  functions  were  generated  to  achieve  convergence  of  the  series  to  at  least  four 
places  of  accuracy  in  both  the  resistive  and  the  reactive  parts  of  the  impedance. 


Fig.  2.  -  Subdivided  oblate  spheroid 


The  problem  was  then  input  to  CHIEF  for  various  subdivision  schemes.  CHIEF  is 
designed  to  take  advantage  of  rotational  and  reflective  symmetry  in  the  geometry,  greatly 
reducing  the  computation  time  in  these  cases.  In  order  to  use  rotational  symmetry,  the 
spheroid  was  first  subdivided  into  strips  resembling  orange  sections.  Each  strip  was 
then  further  subdivided  into  bands.  The  number  of  Gaussian  quadrature  points  for  the 
non- self  Helmholtz  integrals  was  selected  for  both  surface  coordinates.  The  normal 
particle  velocity  was  input  as  unity  for  the  top  half  and  zero  for  the  bottom  half  of  the 
spheroid. 

CHIEF  was  used  to  obtain  the  surface  pre.ssures  Pa  for  each  subdivison.  The  nor¬ 
malized  self  radiation  impedance  coefficient  for  the  top  half  of  the  spheroid,  was 

then  calculated  using 
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^sclf  ”  ^sclf  *^self  “ 


S  Pa 

ggl _ 

n 

pcv 

a  =  l 


(11) 


where  is  the  area  of  the  ath  subdivision,  and  a  r  i  to  n  includes  only  the  subdivisions 
on  the  top  half  of  the  spheroid.  Rseu  and  are  the  resistive  and  reactive  compo¬ 
nents  of  the  impedance.  Here  v  =  1  m/sec  is  the  normal  particle  velocity.  The  nor¬ 
malization  factor 


asjl 


makes  Z^eif  dimensionless. 


Similarly,  the  normalized  mutual  radiation  impedance  coefficient  between  the  two 
halves,  ,  was  calculated  using 


Z  +tV  ^_aan+l 

mutual  "  mutual  ^^'mutual 


N 


(12) 


asn+ 1 


where  a  =  n  +  1  to  N  includes  only  the  subdivisions  on  the  bottom  half  of  the  spheroid. 

The  results  for  this  case  are  summarized  in  Table  1.  The  impedances  calculated 
using  the  harmonic  expansion  are  given  in  the  last  line.  The  column  labeled  "time  in 
reconds"  indicates  the  computation  time  required  to  determine  the  surface  pressures 
using  the  CDC  3800  computer  at  NRL. 

The  first  CHIEF  model  had  12  strips,  6  bands,  used  a  2-point  Gaussian  quadrature 
for  both  surface  coordinates  (denoted  2  x  2),  and  required  36  sec  of  computation  time. 

The  results  are  good  ex-;ept  for  the  mutual  reactance.  As  the  quadrature  is  increased, 
the  accuracy  improves,  with  the  8  x  8  quadrature  giving  surprisingly  good  results  for 
76  sec  of  computer  time. 

If  24  strips  and  12  bands  are  used,  a  2  x  2  quadrature  again  gives  a  poor  value  for 
the  mutual  reactance.  As  the  quadrature  increases,  the  accuracy  again  improves.  In 
general,  the  accuracy  should  be  better  for  increased  subdivision  for  the  same  quadrature. 
However,  possible  random  errors  may  cancel,  yielding  significantly  greater  accuracy 
for  the  coarser  subdivision  scheme.  This  is  apparently  the  case  for  the  model  using  12 
strips,  6  bands,  and  an  8  x  8  quadrature. 

Apparently  good  results  can  be  obtained  for  this  oblate  spheroid  using  rather  crude 
surface  subdivision  and  a  small  number  cf  quadrature  points.  It  is  important  to  note 
that  the  far-field  pressure  pattern  does  not  depend  on  the  radiation  reactance  and  will  be 
extremely  accurate  for  all  of  the  models  except  those  using  a  2  x  2  quadrature.  Because 
of  the  smallness  of  the  mutual  reactance,  the  near-field  pressures  should  also  be  accu¬ 
rate  for  all  but  the  2  x  2  quadrature. 
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Table  1 

Results  for  the  Spheroid  with^  =  0.2,  h  =  1.0 


Number 
of  Strips 

Number 
of  Bands 

Number  of 
Quadrature 
Points 

Time 

(sec) 

*^self 

^Sfl  f 

p 

^mutual 

^mutual 

12 

6 

2x2 

36 

0.2364 

0.5063 

0.1830 

0.06148 

12 

6 

4x4 

49 

0.2371 

0.5254 

0.1760 

0.03996 

12 

6 

8x8 

76 

0.2369 

0.5285 

0.1743 

0.03508 

24 

12 

2x2 

106 

0.2361 

0.5172 

0.1767 

0.04543 

24 

12 

4x4 

1S4 

0.2360 

0.5234 

0.1751 

0.03936 

24 

12 

6x6 

227 

0.2361 

0.5246 

0.1748 

0.03823 

24 

12 

8x8 

338 

0.2361 

0.5249 

0.1748 

0.03791 

Harmonic  expansion 

20 

0.2369 

0.5287 

0.1739 

0.03481 

Example  11.  Thin  Oblate  Spheroid 

Next  consider  the  radiation  from  a  thin  oblate  spheroid  whose  top  half  is  vibrating 
with  unit  velocity  and  whose  bottom  half  is  rigid.  Here  h  =  1.0 ,  and  ^  =  0.02  so  that 
the  ratio  of  major  to  minor  axes  is  very  nearly  equal  to  50.  The  radiation  impedance  of 
the  top  half  of  the  spheroid  and  the  mutual  radiation  impedance  coefficients  between  the 
two  halves  of  the  spheroid  wei'e  calculated  using  a  harmonic  expansion  in  spheroidal 
wave  functions  and  using  CHIEF  with  various  subdivision  schemes.  The  results  are 
given  in  Table  2. 


Table  2 

Results  for  the  Spheroid  with  ^  =  0.02,  h  =  1.0 


Number 
of  Strips 

Number 
of  Bands 

Number  of 
Quadrature 
Points 

— 

^self 

^sclf 

p 

mutual 

mutual 

12 

6 

4x4 

42 

0.2906 

0.3897 

12 

6 

8x8 

65 

0.2513 

0.4247 

0.3160 

24 

12 

4  X  4 

148 

0.2460 

0.4241 

0.3048 

24 

12 

8  X  8 

334 

0.2284 

0.4678 

0.2163 

0.1983 

24 

18 

8x8 

678 

0.2267 

0.4726 

0.2135 

0.1917 

43 

24 

4x4 

770 

0.2235 

0.4714 

Basil 

0.1890 

48 

24 

8x8 

2204 

0.2288 

0.5464 

0.09471 

Harmonic  expansion 

20 

0.2302 

Ugl 

0.06135 

The  results  using  a  model  with  12  strips,  6  bands,  and  a  4  x  4  quadrature  are  very 
poor.  Increasing  the  quadrature  to  8  x  8  does  not  improve  the  results  much.  Increased 
subdivision  gradually  improves  the  results.  However,  even  with  48  strips,  24  bands, 
and  an  8  X  8  quadrature,  considerable  error  occurs  in  the  .autual  reactance.  A  good 
far-fisld  pattern  should  be  obtained  for  this  model,  but  more  subdivisions  and  increased 
computer  time  would  be  required  to  obtain  very  accurate  impedances. 
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Example  III.  Ttvo  Spheres 

To  determine  the  capability  of  CHIEF  regarding  multiple  surfaces,  the  two-sphere 
problem,  whose  geometry  is  shown  in  Fig.  3,  was  considered.  The  spheres  were  chosen 
to  have  ka  =  l.o,  where  a  is  the  radius  of  each  sphere,  and  were  separated  by  a  distance 
equal  to  their  radii.  One  sphere  was  pulsating  uniformly,  while  the  other  one  was  rigid. 
The  two-sphere  problem  has  been  solved  analytically  by  New  (4)  using  expansions  in 
spherical  wave  functions. 


Consider  the  near-field  pressure  magnitude  on  the  axis  of  the  system,  as  shown  in 
Fig,  4.  The  dashed  line  represents  the  l/r  dependence  to  be  expected  if  the  rigid 
sphere  were  not  present.  The  results  of  New  are  represented  by  the  solid  line,  while 
the  dots  represent  the  results  of  CHIEF  using  6  strips,  6  bands,  and  a  4  x  4  quadrature 
on  each  sphere.  The  pressures  are  normalized  to  the  pressure  magnitude  that  would 
exist  on  the  surface  of  the  pulsating  sphere  if  the  rigid  sphere  were  not  present.  The 
maximum  error  was  less  than  6%.  When  10  strips,  10  bands,  and  a  4  x  4  quadrature 
were  used,  the  maximum  error  was  reduced  to  less  than  1%. 

The  far-field  pressure  pattern  in  the  plane  bisecting  the  spheres  was  also  calculated 
using  the  CHIEF  model  with  6  strips,  6  bands,  and  a  4  x  4  quadrature  on  each  sphere. 
Figure  5  gives  the  pattern  as  a  function  of  the  polar  angle  o.  Again  the  solid  line  repre¬ 
sents  the  results  of  New,  while  the  CHIEF  results  are  represented  by  dots.  As  expected, 
considering  the  good  accuracy  of  the  near-field  results,  the  far-fitld  is  extremely 
accurate.  The  total  computation  time  for  the  two-sphere  problem  using  the  coarser 
subdivision  scheme  was  126  sec. 


GEOMETRICAL  OPTIONS 

CHIEF  is  not  restricted  to  a  specific  coordinate  geometry.  Any  convenient  coor¬ 
dinate  representation  may  be  used  to  deocribe  the  surface.  However,  the  free-space 
Green's  function  and  its  normal  deru^ative  are  described  internally  in  CHIEF  in  terms 
0''  a  Cartesian  coordinate  system.  Therefore,  input  subroutines  are  required  which  con¬ 
tain  conversion  formulas  giving  the  Cartesian  coordinates  of  an  arbitrary  surfaxe  point 
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Fig.  4  -  Near-field  pressure  distribution  on  the 
_xls  of  two  spheres,  one  pulsating  and  one  rigid 


Fig.  5  -  The  far-field  pressure  distribution  for 
two  spheres,  one  pulsating  and  one  rigid 

as  well  as  the  Cartesian  components  of  the  normal  vector  to  the  surface  at  that  point. 
The  area  element  associated  with  the  surface  coordinates  is  also  required.  The  entire 
closed  surface  is  separated  into  regions  such  that  each  region  is  describable  in  terms  of 
a  single  geometry.  Each  region  is  assigned  an  integer  index  that  corresponds  to  the 
appropriate  conversion  formulas  in  the  input  subroutines. 

In  an  effort  to  increase  the  utility  of  CHIEF,  conversion  formulas  have  been  added 
to  the  input  subroutines  so  that  a  wide  rai^e  of  geometrical  options  is  available.  These 
options,  which  are  listed  in  Table  3,  are  discussed  below.  Included  are  formulas  giving 
the  Cartesian  coordinates  (x,  y,  z),  the  Cartesian  components  of  the  imit  normal  Cn, 
and  the  magnitude  of  the  area  element  dS  for  any  surface  point  (u,  v)  (Ref.  5).  The 
notation  (u.  v)  a  (n,  b)  means  that  a  and  b  are  the  two  surface  coordinates.  Any  of 
the  surfaces  can  be  translated  or  rotated  by  modifying  the  formulas  giving  the  Cartesian 
coordinates.  For  example,  option  15  describes  the  outside  surface  of  a  sphere  centered 
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at  (x,  y,  z)  =  (0,  0,  0).  To  describe  a  sphere  centered  at  (x,  y,  z)  =  (c,  o.  O),  change 
the  equation  for  x  in  terms  of  the  surface  coordinates  (u,  v)  a  (0,  4>)  to 
X  =CC(33)  sinu  cos  v  +  c. 

i.  yz  plane  with  Cartesian  coordinates,  (u,  v)  a  (y,  z),  normal  in  ^  x  direction, 

X  =  Constant  =  CC(1) 
y  =  u 
z  =  V 

=  Sx 

dS  =  du  dv 

Si.  yz  plane  with  Cartesian  coordinates,  (u,  v)  s  (y,  z),  normal  in  -x  direction. 

X  =  Ctxistant  =  CC(2) 

y  =  u 

z  =  V 

0  =  -6 
n  X 

dS  =  du  dv 

3.  xz  plane  with  Cartesian  coordinates,  (u,  v)  s  (x,  z),  normal  in +  y  direction. 

X  =  u 

y  =  Cbnsiant  =  CC(3) 


Table  3 

Geometrical  Options  Added  to  CHIEF 


Surface 

Direction  of 
Normal  to 
Surface 

Surface 

Direction  of 
Normal  to 
Surface 

1.  yz  plane 

.X 

11.  Outside  of  circular 

Outward 

2.  yz  plane 

-X 

cylinder 

3.  xz  plane 

^•y 

12.  Inside  of  circular 

Inward 

4.  xz  plane 

-y 

cylinder 

5.  xy  plane 

1*  Z 

13,  Outside  of  elliptical 

Outward 

6.  xy  plane 

—  z 

cylinder 

7.  xy  plane,  polar 

+  7. 

14.  Inside  of  elliptical 

Inward 

coordinates 

cylinder 

8.  xy  plane,  polar 

-Z 

15.  Sphere 

Outwai  d 

coordinates 

16.  Oblate  spheroid 

Outward 

9.  xy  plane,  elliptical 

^  z 

17.  Prolate  spheroid 

Outward 

coordinates 

18.  Toroid 

Outward 

10.  xy  plane,  elliptical 

coordinates 

A.L.  VANBUREN 


e  =  e 
n  y 


(Is  =  du  dv 


4.  xz  plane  with  Cartesian  coordinates,  (u,  v)  a  (x.  z),  normal  in  -y  direction . 


y  =  Constant  =  CC(4) 


dS  =  du  dv 


5.  xy  plane  with  Cartesian  coordinates,  (u,  v)  a  (x,  y),  normal  in  +  z  direction. 


z  =  Constant  =  CC(S) 


dS  =  du  dv 


6.  xy  plane  with  Cartesian  coordinates,  (u,  v)  a  (x,  y) ,  normal  in  -z  direction. 


z  =  Constant  =  CC(6) 


6  =  -e 

11  z 


dS  =  du  dv 


7,  xy  plane  with  polar  coordinates  (r.  O),  0^  r  <«»,  0  <  277,  (u,  v)  s  (r.  0) ,  normal 

in  +z  direction. 

In  the  xy  plane  the  curve  of  constant  r  is  a  circle  of  radius  r  centered  at  (X.  y) 
=  (0,  0),  and  the  curve  of  constant  is  a  semi-infinite  line  originating  at  (x,  y) 
=  (0,  0). 

X  =  U  (X5S  V 

y  =  u  sin  v 
z  =  Const.mt  =  CC(7) 


(IS  =  u  du  dv 


NRL  REPORT  7160 


11 


xy  plane  with  polar  coordinates  (r,  0),  r  <  ®,  Qi  0  <  2v,  (vi,  v)  a  (r,  0) 
normal  in-z  direction. 

In  the  xy  plane  the  curve  of  constant  r  is  a  circle  of  radius  r  centered  at 
(x,  y)  =  (0,  0),  and  the  curve  of  constant  5  is  a  semi-infinite  line  originating  at 
(X,  y)  =  (0,  0)  • 


X  =  u  cos  V 
y  =  u  sin  v 
2  =  Constant  =  CC(8) 


dS  =  u  du  dv 


xy  plane  with  elliptical  coordinates  (m.  1  0  <  «/-  <  2v,  (u,  v)  a  {fx,  i/i), 

normal  in  +z  direction. 

In  the  xy  plane  the  curve  of  constant  /n  is  an  ellipse  of  interfocal  distance 
2  CC(9)  centered  at  (x,  y)  =  (0,  O),  and  the  curve  of  constant  is  a.  hyperbola  which 
is  orthogonal  to  the  family  of  ellipses  for  the  same  value  of  cc(9). 


X  =  CC(9)  u  cos  V 
y  =  CC(9)  (u^  -1)^  sin  V 
2  =  Constant  =  CC(IO) 
e  =  8 

z 

dS  =  CC(9)^  (u^  -  1)"^  (u^  -cos^  v)  du  dv 


xy  plane  with  elliptical  coordinates  (m.  i/”  2-n,  (u,  v)  a  {fx,  <//), 

normal  in  -z  direction. 

In  the  xy  plane  the  curve  of  constant  m  is  an  ellipse  of  interfocal  distance 
2  CC(  12)  centered  at  (x,  y)  =  (O,  O),  and  the  curve  of  constant  >/'is  a  hyperbola 
which  is  orthogonal  to  the  family  of  ellipses  for  the  same  value  of  CC(12) . 


X  =  CC(12)  u  cos  V 
y  =  CC(12)  (u^ -ly-^  sin  V 
2  =  Constant  =  CC(13) 


dS=  CC(12)^  (u^  -  1)~^  (u^-cos^v)  du  dv 
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11.  Outside  surface  of  a  circular  cylinder  along  the  z  axis,  (u,  v)  2  (z,  0)  -w  <  z  <  a-, 
0  <0  <  2v,  normal  outward  from  the  cylinder.  The  radius  of  the  cylinder  is  a  con¬ 
stant  given  by  cc(l5) . 

On  the  cylindrical  surface  the  curve  of  constant  z  is  a  circle  of  radius  CC(15) 
centered  at  (x,  y,  z)  =  (O,  0.  2),  and  the  curve  of  constant  t?  is  a  line  parallel  to 
the  z  axis. 


X  =  CC(15)  cos  V 
y  =  CC(15)  sin  v 
z  =  u 

=  Cjj  cos  V  +  sin  v 
dS=  CC(IS)  tki  dv 

12.  Inside  surface  of  a  circular  cylinder  along  the  z  axis,  (u,  v)a  (z,  6),  -m<  z  <  <0,0  <e  <  2n, 
normal  inward.  The  radius  of  the  cylinder  is  a  constant  given  by  CC(lb). 

On  the  cyiindrical  surface  the  curve  of  constant  z  is  a  circle  of  radius  CC(  16) 
centered  at  (x,  y,  2)  =  (O,  0.  z),  and  the  curve  of  constant  0  is  a  line  parallel  to  the 
z  axis. 


X  =  CC(16)  cos  V 
y  =  CC(16)  sin  v 

Z  =  VI 

“  “^x  ''  ““y  *'■ 

dS  =  CC(16)  di  dv 

13.  Outside  surface  of  an  elliptic  cylinder  along  the  z  axis,(u.  v)  3  (z,  1/1),  -  co  <  z  <  co, 
0  <  i/(  <  27r,  normal  outward  from  the  cylinder.  The  interfocal  distance  of  the  ellip¬ 
tic  cross  section  is  given  by  2CC(17).  The  surface  is  one  of  constant  /i  =  CC(18), 
where  1  <  /z  <  <». 

On  the  cylindrical  surface  the  curve  of  constant  z  is  an  eliipse  of  interfocai 
distance  2CC(17)  centered  at  (x.  y,  z)  =(0,  0,  z),  and  the  curve  of  constant  i/>  is  a 
line  parallel  to  the  z  axis,  one  of  the  two  lines  produced  by  the  intersection  of  a 
hyperbolic  cyiinder  and  the  eiiiptic  cylinder. 

X 

y 

z 


dS 


=  CC(17)  CC(18)  cos  V 
=  CC(17)  (CC(18)2  -  1]«  sin  v 
=  11 

I  ^ 


cccib)*^-  1 


CC(18)2-cos2 


V 


cos  V  +  c.. 


CC(18) 


ICC(18)2-  cos^  v]’ 


=  CC(17)  ICC(18)2  -cos^  v]  du  dv 
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14.  Inside  surface  of  an  elliptic  cylinder  along  the  z  axis,  (u,  v)  a  (z, ./)),-  ®<  z  <  a>, 

0  <  ■/-  <  277 ,  normal  inwai'd.  The  interfocal  distance  of  the  elliptic  cross  section  i.s 
given  by  2CC(25).  The  surface  is  one  of  constant  /i  =  cc<'26),  where  i  <  (i<  <j>. 

On  the  cylindrical  surface  the  curve  of  constant  z  is  an  ellipse  of  interfocal 
distance  2CC(25)  centered  at(x,  y.  z)  =  (O.  o.  z),  and  the  curve  of  constant  ■/-  is  a  line 
parallel  to  the  z  axis,  one  of  the  two  lines  produced  by  the  intersection  oi  a  hyper¬ 
bolic  cylinder  and  the  elliptic  cylinder, 

X  =  CC(25)  CC(26)  cos  v 

y  =  CC(25)  [CC(26)2  -  t]  **  sin  v 


g  r  CC(26)^-1 
*  [cC(26)2-  cos^vj 


.  C'C(  26) 

COS  V  -  ,  —  , — i — -  -  ■  sin  v 

ICC(26)2_cos2  v]’^ 


dS  =  CC(25)  [CC(26)2  -cos^v]  du  dv 


15.  Outside  surface  of  a  sphere,  (u.  v)  k  (0.  <(>),  Q  <d  <v,  O  <(j><  2ti,  normal  outward 
from  the  sphere.  The  radius  of  the  sphere  is  a  constant  given  by  CC(33). 

On  the  spherical  surface  the  curve  of  constant  ^  is  a  circle  of  radius 
cc(33)sin0  centered  at  (x,  y,  t)  =  (O,  o.  CC(33)  cos^),  and  the  curve  of  constant  <i> 
is  a  half  circle  between  the  poles  of  the  sphere. 


X  =  CC(33)  sin  u  cos  v 
y  =  CC(33)  sin  u  sin  v 
z  =  CC(33)  cos  u 

=  Sj,  sin  u  cos  V  +  Cy  sin  u  sin  v  +  6^  cos  u 
dS  =  CC(33)  sin  u  du  dv 

16.  Outside  surface  of  an  oblate  spheroid,  (u,  v)  s  (t?,  0),  0  <  t?  <  77,  0  <<t><  2v,  normal 
outward  from  the  spheroid.  The  interfocal  length  of  the  spheroid  is  a  constant 
given  by  2CC(34).  The  surface  is  one  of  constant  ^  =  CC(35),  where  0  <  ^  <  ®. 

On  the  spheroidal  surface  the  curve  of  constant  is  a  half  ellipse  containing 
the  poles  of  the  spheroid,  and  the  curve  of  constant  77  is  a  circle  of  radius 
CC(34)  [^2  +  1]**  sin  77  centered  at  (x,  y,  z)  =  (0,  0,  CC(34)  cx:(35)  cos  77).  The  Z  axis  is 
the  axis  of  symmetry. 


X  =  CC(34)  tCC(3S)^  +  1]**  sin  11  cos  v 
y  =  CC(34)  (CC(35)^  +  I]'*  sin  u  sin  v 
z  =  CC(34)  CC(35)  cos  u 
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^C(3S)  tCC(35)^  >  cos^  u]  **  sin  u  cos  v 
+  ‘■‘y  CC(3S)  [CC(3S)^  +  cos^  uj  ^  sin  u  sin  v 
+  Cj  [CC(3S)^  +  1]  **  iCC(35)“  +  cos^  u]~**  cos  u 
~  CC(34)^  (CC(3S)^  +1]^  lCC(35)^  +  cos^  u]**  sin  u  du  dv 

17.  Outside  surface  of  a  prolate  spheroid,  (u,  v)  a  (tj,  i^),  o  <-n  <n,  o  <q}<  2v ,  normal 
outward  from  the  spheroid.  The  interfocal  length  of  the  spheroid  is  a  constant 
given  by  2CC(42) .  The  surface  is  one  of  constant  ^  =  CC(43),  where  i  «  ^  <  oo . 

On  the  spheroidal  surface  the  curve  of  constant  is  a  half  ellipse  containing 
the  poles  of  the  spheroid,  and  the  curve  of  constant  -n  i.,  a  circle  of  radius 
CC(42)  [4^  -  1]*^  sin  t)  centered  at  (x,  v.  2)  =  (O.  0,  CC(42)  CC(43)  cos  7)),  The  z 
axis  is  the  axis  of  symmetry . 

X  =  CC(42)  [CC(43)2  -  1]  **  sin  u  cos  v 
y  =  CC(42)  [CC(43)2-1]“  sin  u  sin  v 
2  =  CC(42)  CC(43)  cos  u 

CC(43)  (CC(43)^  -  cos^  u]~**  sin  u  cos  v 
+  Sy  CC(43)  [CC(43)^  —  co3^  u)”**  sin  u  sin  v 
+  S,  [CC(43)2  -  1]^  [CC(43)2  -  cos*  u)"**  cos  u 
dS=  CC(42)*  [CCf43)* -1]'*  (CC(43)* -cos*  u3^  sin  u  du  dv 

18.  Outside  surface  of  a  toroid,  (u,  v)  a  (tj,  <^),  -n  <  ri  <  n,  0  <  <f>  <  2tt,  normal  outward 
from  the  toroid.  The  toroid  is  characterized  by  two  radii  cc(49)  and  CC(50).  The 
ratio  of  these  two  radii  CC(5l),  1  <  CC(51)  <  ® ,  forms  an  orthogonal  system  with  tj 
and  (j}  and  is  a  constant  over  the  toroidal  surface.  To  completely  characterize  the 
surface,  a  second  constant  CC(52)  =  (CC(SO)* -00(49)*]**  is  defined. 

Thus  the  surface  may  be  defined  by  giving  either  the  constants  00(49)  and 
00(50)  or  the  constants  00(51)  and  00(52).  Because  the  two  radii  are  easier  to 
visualize,  00(49)  and  00(50)  will  be  the  input  for  CHIEF. 

On  the  toroidal  surface  the  curve  of  constant  is  a  circle  produced  by  the  inter 
section  of  a  half  plane  containing  the  z  axis  and  the  toroid,  and  the  curve  of  constan 
TJ  is  a  circle  produced  by  the  intersection  of  a  spherical  bowl  and  the  toroid.  The 
2  axis  is  the  symmetry  axis  for  the  toroid. 


X  =  (00(50)* -00(49)*]  cos  v/(00( 50) -00(49)  cos  u] 
y  =  (00(50)* -00(49)*]  sin  v/(00(50) -00(49)  cos  u] 

2  =  00(49)  (00(50)* -CO (49 )*]■'*  sin  ii/(00(50) -00(49)  cos  u) 
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6„  =  -  [CC<49)  -CC(SO)  cos  u]  cos  v/[CC(S0)  -  CC(49)  cos  u] 

-  gy  (CC(49)  -CC(SO)  cos  u]  sin  v/tCC(S0)  -CC(49)  cos  u] 
+  g^  [CC(S0)2  -00(40)2]*^  sin  u/(CC(S0)  -  CC(49)  cos  u] 
dS  =  CC(49)  [CC(50)2  -CC(49)^]'V[CC(S0)  -  CC(49)  cos  u)^  du  dv 


The  input  subroutines  providing  for  these  geometrical  options  are  listed  in  the 
appendix.  Note  that  the  two  subroutines  CCOORD  and  CCUNMD  already  exist  in  CHIEF, 
CCOORD  (U,  V,  NEQN,  Al,  A2,  A3)  is  called  in  CHIEF  whenever  only  the  Cartesian 
coordinates  (Al,  A2,  A3)  are  desired  for  a  surface  point  (U,  V).  The  index  NEQN  refers 
to  the  particular  geometrical  option.  The  subroutine  CCUNMi)  (U,  V,  NEQN,  AXl, 

AX2,  AX3,  ANl,  AN2,  AN3,  RMAGD)  is  called  whenever  all  of  the  surface  information 
is  desired;  i.e,,  it  provides  the  Cartesian  coordinates  (AXl,  AX2,  AX3),  the  Cartesian 
components  of  the  unit  normal  to  the  surface  (ANl,  AN2,  AN3),  and  the  magnitude  of 
the  surface  area  element  RMAGD. 

Note  that  the  set  of  constants  Cl,  C2,  C3, , . ,  in  CHIEF  have  been  replaced  by  an 
array  of  constants  CC(IOO).  The  declaration  COMMON/ALLC/CC(100)  must  be  sub¬ 
stituted  in  the  main  program  and  in  CCOORD  and  CCUNMD  for  the  declaration 
COMMON /ALLC /Cl,  C2,  C3, . . .  These  constants  are  quantities  that  do  not  vary  over 
the  surface  and  must  be  input  or  calculated  in  the  main  program.  Every  combination 
of  constants  appearing  in  the  conversion  formulas  is  calculated  and  stored  in  CC(IOO). 
This  reduces  the  computation  time  considerably,  since  the  calculation  is  only  performed 
once  instead  of  every  time  the  subroutines  CCOORD  and  CCUNMD  are  called. 

Most  of  the  input  constants  refer  to  a  distance  and  are  described  in  terms  of  a  unit 
of  length  called  WAVE=  lA  =  k/lv  =  c/w  ,  However,  some  of  the  input  constants  in 
options  13,  14,  16,  and  17  are  pure  numbers  and  represent  ratios  of  distances.  E'.efer- 
ence  to  the  previous  asscription  of  the  geometrical  options  will  indicate  which  parameters 
the  input  constants  represent. 

Option  19  in  subroutine  CCOORD  provides  interior  points  and  should  be  modified  to 
fit  the  specific  geometry  whenever  interior  points  are  required. 


SUMMARY 

The  following  observations  may  aid  the  potential  user  of  CHIEF. 

1.  CHIEF  evaluates  the  solution  to  a  boundary-value  problem.  This  will  be  a  realistic 
transducer  model  when  the  velocities  are  known;  for  example,  when  velocity  control 
exists , 

2.  CHIEF  is  a  very  flexible  program,  allowing  the  user  a  wide  range  of  options. 
Reducing  it  to  a  production  program,  capable  of  being  used  by  casual  acquaintances, 
would  remove  this  flexibility.  If  only  a  specific  geometry  is  of  interest,  the  user  is 
advised  to  write  a  CHIEF-lilse  program  taking  advantage  of  the  peculiarities  of  that 
geometry.  This  will  l)e  more  economical  in  the  long  run,  if  the  program  is  to  be  used 
extensively.  If  the  frequencies  of  interest  are  w'ell  below  the  lowest  critical  frequency, 
a  simple  source  method,  where  the  pressure  is  obtainable  in  terms  of  a  single  integral, 
might  be  considered . 

3.  If  rotational  or  reflective  symmetry  is  not  present,  the  computation  times  increase 
significantly.  The  possibility  of  using  rotational  or  reflective  symmetry  should  be  ex¬ 
amined  carefully.  If  both  are  applicable,  it  is  usually  best  to  use  rotational  symmetry. 
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4.  The  results  obtained  using  CHIEF  approach  the  correct  solution  as  the  subdivision  is 
increased.  Surprisingly  gcod  results  can  be  obtained,  especially  for  the  far  field,  using 
a  relatively  crude  subdivision  scheme.  However,  objects  that  are  thin  may  require  an 
unusually. fine  subdivision.  Also,  as  the  acoussdc  size  is  increased,  the  number  of  sub¬ 
divisions  reqiiired  to  give  accurate  results  will  also  increase.  However,  this  limitation, 
which  is  fundamental  to  the  finite-element  method,  will  be  less  restrictive  in  the  future 
as  computers  are  improved. 

5.  This  report  describes  new  input  subroutines  to  CHIEF  that  provide  for  a  wide  range 
of  geometrical  options.  Most  of  the  standard  geometries  have  been  included.  Surface.'? 
of  complex  objects  such  as  a  free-flooded  ring  can  be  described  by  combining  two  or 
more  of  the  options.  The  user  can  easily  add  additional  geometries  tliat  are  desired. 
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t4  Al=CC(a<?)»COS(V) 
A2=CC(30)*SIN{V) 

A3=U 

RETURN 

15  RSINU=CC(33)*SIN(U> 

A1 =RSINU»COS( V) 
A2=RSlNU*StN< V) 
A3=CC(33)*C0S(U) 

RETURN 

16  CCSINU=CC(39)*S1M(U) 
A1=CCS1NU»C0S<V) 
A2=CCSINU*SIN< V) 
A3=CC(37)*C0S(U> 

return 

17  CCSINU=CC(A7)*S1N<U) 
A1=CCSINU*C0S<V) 
A2=CCSINU*S1N< V) 
A3=CC<45)*C0S(U) 

RETURN 

18  FAC3CC<54)/'CC<51 J-COS(U) ) 
A1 =FAC*COS( VI 
A2=FAC*SIN(V1 
A3=FAC*SIN(U)*CC<S5) 

RETURN 

19  A1=0. 

A2=U 

A3=CC(60l 

RETURN 

END 


SUOROUTINE  CCUNMD<U«V«NEQN,AX1 .AX2*AX3«ANl ♦ AN2« AN3*RMAGD> 
COMMON/ALLC/CC< 100 ) 

GOTO  ( 1 «2«3«4.5»6«7.8»9* 10« 1 1 ♦ 12* 13* 14. I5t 16* 17* 18>  NEON 

1  AN1=1. 

AN2=0, 

AN3=0, 

RMAGOr 1 . 

AX1=CC( 1 ) 

AX2=U 

AX3=V 

return 

2  AN1=-1 . 

AN2=0, 

AN3=0. 

RMAGn=l • 

AXinCCl?) 

AX2=U 

AX3=V 

RETURN 

3  AN1=0. 

AN2=1 . 

AN3=0, 

RMAGO- 1 • 

AX|  =U 
AX2=CC(3) 

AX3=V 

RETURN 

4  AN1=0, 

AN?  =  -1  . 

AN3=0. 

RMAGr)=  1  • 

AX  1  =U 
AX2=CC(4) 
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AX3=V 

return 

*5  AN1=0, 

AN2=0, 

AN3= ! . 

.PM  A  fin  =  1  • 

AX1=U 

AX2=V 

AX3=CC(5) 

RETURN 

6  ANlrO, 

AN2=0. 

AN3=-1 . 

RMAGn=l 4 
AX1=U 
AX2=V 
AX3=CC(6) 

RETURN 

7  AN1=0, 

AN2=0, 

AN3=1 4 

RMAGn=U 

AXl=U*COfi(V) 

AX2=U#S!N(V) 

AX3=r,C(7) 

RETURN 

a  AN1=04 
AN2=04 
AN3=-1 4 
RMAfiRsU 
AXl=U*COS(V> 

AX2sU#StN{V) 

AX3=CC<8) 

RETURN 

<?  COSV=COS<V) 

COSVS=COSV»COSV 

UU=U»U 

SO=SORT(UU-l 4 ) 

AN1=04 

AN2=04 

AN3=l4 

RMAGD=CC( I J )»(UU-COSVS)/SQ 
AX 1 =CC ( 9 ) *U*COSV 
AX2=CC  <  9 ) »SQ*SORT ( 1 4 -COSVS ) 
AX3=CC( 10) 

RETURN 

10  COSV=COS(V) 

UU=U*U 

cosvs=cosv»cosv 

SQ=S0RT<UU-1 4  > 

AN1=04 
AN2=0  4 
AN3=-1 4 

RMAGn=CC { 1 A ) * « UU-COSVS ) /SO 
AX1=CC{ 12)#U*C0SV 
AX2=CC ( 1 2 ) #SO*SORT ( 1 4 -COSVS ) 
AX3=CC( 13) 

RETURN 

11  AN1=C0S(V) 

AN2=S1N<V) 

AN3=04 

RMAGr)=CC(  16) 

AX1=CC( 15)*AN1 
AX2=CC( 15)»AN2 
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AX3=U 

RETURN 

12  ANI=-COS(V) 

AN2=-S1N( V) 

ANSsO, 

RMAGD=CC( 16) 

AX1=-CC( 161*An1 
AX2=-CC( 16)»AN2 
AX3=U 
RETURN 

13  COSV=COS(V) 

COSVS=COSV*COSV 
SINV=S0RT( 1,-COSVS) 

OENS=l .n/SORT<CC< 19)-COSVS) 
AN  1  =CC  {  2C 1  »r>ENS*COSV 
AN2=cc( ia)*nENs*s:Nv 

AN3=0. 

RMAGD=CC( 17) /DENS 

AXlsCC(21 )*COSV 

AX2=CC<2P)*AtNV 

AX3=U 

RETURN 

lA  COSV=COS<V) 

COSVS=COSV»COSV 
S 1 NV -  SORT ( 1 . -COSVS ) 

OENS= 1 . /SORT ( CC ( 27 ) -C03VS ) 
AN 1 =-CC ( 2fi ) *OFNa»COSV 
ANP=-r.C  (  26 )  *OFNS*Sl  NV 
AN3=0. 

RMAGD=CC ( 25 ) /DENS 

AX1=CC(29)»C0SV 

AX2=CC(30)»StNV 

AX3sU 

RETURN 

15  SINU=SIN(U) 

ANl=SINU»COS(V) 

AN2sSlNU*SIN(V) 

AN3=COS('J) 

RMAGD=CC(33)*SINU 
AXl rCCla3)»ANl 
AX2=CC<33)*AN2 
AX3sCC<33)*AN3 

return 

16  COSV=COS<V) 

SlNV=StN(V) 

COSU=COS(U) 

cosur.=cosu»cosu 

SINU=SQRT{ 1 .-COSUS) 

FAC=1 ./SORT<CC<36)+COSUS) 

CCSINU=CC(39) *SINU 

GAC=FAC*CC(35)»SINU 

ANl =GAC*COSV 

AN2=GAC*SINV 

AN3=FAC»CC ( 36 ) »COSU 

RMAGD=CC ( AO ) »Sl NU/FAC 

AXl =CCSINU»COSV 

AX2=CCSINU*S1NV 

AX3=CC<37)*C0SU 

RETURN 

17  COSV=COS(V) 

S1NV=SIN(V) 

COSU=COS<U) 

cosus=cosu*cosu 


I 

I' 
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SINU=S0RT( 1 ,-COSUS) 

FAC= 1  * /SORT ( CC ( 4 A ) -COSUS ) 

CCS INU=CC ( 47 ) »S I NU 

GAC=FAC*CC(4S)*SJNU 

ANI=GAC;*C0SV 

AN?=GAC»SIMV 

AN3=F AC*CC { 46 » *COSU 

RMAGD=CC<48>*SINU/FAC 

AXl=CCSINU»COSV 

AX2=CCSINU»SINV 

AX3=CC(45)»C0SU 

return 

IR  COSV=COS<V) 

SJNV=SIN( V) 

COSU=COS(U) 

StNU=SIN(U) 

BAC- 1 . / ( CC ( 5 1 ) -casu ) 

CAC=BAC*< 1 .-CC*^! )*COSU> 

FAC=nAC»CC(<54  ) 

ANl=-CAC1‘COSV 

AN2=-CAC*S!NV 

AN3=  CC(53)»SINU»RAC 

RMAGD=CC ( 56 ) *BAC*SAC 

AXl=FAC*COSV 

AX2=FAC#SJNV 

AX3=FAC»S !NU»CC ( 5B ) 

RETURN 

END 


1  CC< J )=1 .*WAVE 

2  CC(2)=l4*WAVE 


3 

CC 

(3) 

=1 «*WAVE 

4 

CC 

<4) 

=1«*WAVE 

5 

CC 

(5) 

=1.»WAVE 

6 

CC 

(6) 

=1 .*WAVE 

7 

CC 

<71 

=»57S»WAVE 

8 

CC 

(8) 

=-.575*WAVP 

9 

CC 

(91 

=1 .»wAvr 

CC 

(  10 

1  =  1  . 

CC 

(  1  1 

)=CC(9)*CC(9) 

10 

CC 

<  12 

1  -  : .#WAVE 

CC 

(  13 

1  =  1 . 

CC 

<  14 

)=CC( 12)*CC(  1 

2) 

1  1 

CC 

(15 

1  =  1 .00*WAVE 

12 

CC 

(  16 

) = 1 .85*WAVe 

13 

CC 

(  17 

>=1 .*WAVE 

CC 

(  18 

1  =  1  . 

CC 

(  19 

)  =  cc{  ianfcc(  1 

8) 

CC 

(20 

)=SORT(CC( 19) 

-1  .  ) 

CC 

(21 

)=CC( 17>*CC( 1 

8) 

CC 

(22 

)=CC(  17)  ftCC(20) 

22 
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14  ccisS) =i .*wAve 
CC(36)=1 . 

CC(27)=CC(26)»CC(26) 

CC<28)=SQRT(CC(?7)-1 .) 
CC(29)=CC<28)*CC(26» 

CC ( 30 ) =CC  <  25 ) «CC  <  28 1 

15  CC<33)=1 .*WAVF 

16  CC(341=l.*WAVe 
CC ( 35 1 = . 2 

CC(36)=CC(35i*CCt35) 

CC<37l=CC(34)*CC<35) 

CC (38) =S0RT<  CC( 36) +1 » ) 
CC(39)=CC<34)*CC(38) 

CC  ( 40  )  =CC.  ( 39 )  *CC  (  34  ) 

17  CC(42)=1 .♦WAVE 
CC(43)=1 . 

CC(44)=CC(43)»CC(43) 

CC(4S)=CC(42)»CC(43) 

CC(46)=S0PT(CC<44)-1.) 

CC(47)=CC(42)*CC<46) 

CC(48)=CC<47)*CC(42) 

18  CC(49)=» la^WAVE 
CC<50)=1.86»WAVE 
CC(51 1=CC(50) /CC(49) 

CC ( 52 ) s SORT ( CC ( SO ) ♦CC  <  50 ) -CC  <  49  > ♦CC ( 49 ) ) 
CCi53)=CC(52)/CC<49) 

CC<54)=CC(52)«CC(S3) 

CC<55)=1 ./CC(53S 
CC<56)»CC(54)^CC<52) 
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